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(N '. Abstract 

Quantum molecular dynamic simulations have been employed to study the equation of state 
(EOS) of fluid helium under shock compressions. The principal Hugoniot is determined from EOS, 
where corrections from atomic ionization are added onto the calculated data. Our simulation results 
indicate that principal Hugoniot shows good agreement with gas gun and laser driven experiments, 
and maximum compression ratio of 5.16 is reached at 106 GPa. 
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I. INTRODUCTION 



High-pressure introduced response of materials, which requires accurate understandings 
of the thermophysical properties into new and complex region, has gained much scientific 
interest recently [l]. The relative high temperature and high density are usually referred as 
the so-called "warm dense matter" (WDM)- a strongly correlated state, where simultaneous 
dissociations, ionizations, and degenerations make modelling of the dynamical, electrical, 
and optical properties of WDM extremely challenging WMD, which provides an ac- 
tive research platform by combining the traditional plasma physics and condensed matter 
physics, usually appears in shock or laser heated targets inertial confinement fusion 4], 
and giant planetary interiors js|. 

Next to hydrogen, helium is the most abundant element in the universe, and physical 
properties of warm dense helium, especially the EOS, are critical for astrophysics [6|, 7 . 
For instance, the structure and evolution of stars, White Dwarfs, and Giant Planets 8- 



llj], and therefore the understanding of their formation, depends sensitively on the EOS 



of hydrogen and helium at several megabar regime. For all planetary models, accurate 
EOS data are essential in solving the hydrostatic equation. As a consequence, a series of 
experimental measurements and theoretical approaches have been applied to investigate the 
EOS of helium. Liquid helium was firstly single shocked to 15.6 GPa using two stage light 
gas gun by Nellis et al, then double shocked to 56 GPa, and the calculated temperature are 
12000 and 22000 K respectively 12|. Maximum compression ratio (f] ma x~Q) was achieved by 
laser driven shock experiments with the crossover pressure around 100 GPa 13j . However, 
soon after that, Knudson et al. have modified rj max to be 5.1 Theoretically, since the 
ionization equilibrium is not interfered with the dissociation equilibrium between molecules 
and atoms, helium, which is characterized by monoatomic molecule and close shell electronic 
structure, is particularly suitable for the investigation of the high pressure behavior under 
extreme conditions. Free energy based chemical models by Ross et al. 151, Chen et al. 

n □ 

|16| . and Kowalski et al. [17J have been used to investigate the principal Hugoniot of liquid 
helium, and the results are accordant with gas gun experiments. However, considerable 
controversies have been raised at megabar pressure regime, especially since the data was 
probed by laser shock wave experiments. Interatomic potential method predict r] max =4: 
for shock compressed helium 15], whereas, the EOS used by the astrophysical community 
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from Saumon et al. (SCVH) [8], path-integral Monte Carlo (PIMC) 18], and activity 
expansion (ACTEX) 19[ calculations provide an increase in compressibility at the beginning 
of ionization. For the initial density of po = 0-1233 g/cm 3 , SCVH and ACTEX simulation 
results indicate the maximum compression ratio lies around 6 at 300 GPa and 100 GPa 
respectively {§], Q, while PIMC calculations suggest r] max =5.3 near 360 GPa [l^. 

On the other hand, quantum molecular dynamic (QMD) simulations, where quantum 
effects are considered by the combinations of classical molecular dynamics for the ions and 



density functional theory (DFT) for electrons, have already been proved to 



}e successful in 



describing thermophysical properties of materials at complex conditions |20l . |21| . However, 
the DFT based molecular dynamic simulations (with or without accounting for excited 
electrons) do not provide reasonable results at the ionization region, mainly because the 
atomic ionization is not well defined in the framework of DFT. Considering these facts 
mentioned above, thus, in the present work, we applied the corrected QMD simulations 
to study shock compressed helium, and the EOS, which is compared with experimental 
measurements and different theoretical models, are determined for a wide range of densities 
and temperatures. The calculated compression ratio is substantially increased according to 
the ionization of atoms in the warm dense fluid. 



II. COMPUTATIONAL METHOD 



The Vienna Ab-initio Simulation Package (VASP) 22|, |23J, which was developed at the 
Technical University of Vienna, has been employed to perform simulations for helium. The 
elements of our calculations consist of a series of volume-fixed supercells including N atoms, 
which are repeated periodically throughout the space. By involving Born-Oppenheimer 
approximation, electrons are quantum mechanically treated through plane-wave, finite- 
temperature (FT) DFT 24|, where the electronic states are populated according to Fermi- 
Dirac distributions at temperature T e . The exchange-correlation functional is determined 
)y generalized gradient approximation (GGA) with the parametrization of Perdew-Wang 91 



251 ] . The ion-electron interactions are represented by a projector augmented wave (PAW) 



pseudopotential [2a]. Isokinetic ensemble (NVT) is adopted in present simulations, where 
the ionic temperature Tj is controlled by Nose thermostat 27|, and the system is kept in 
local equilibrium by setting the electron (T e ) and ion (Tj) temperatures to be equal. 
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The plane-wave cutoff energy is selected to be 700.0 eV so that the pressure is converged 
within 3% accuracy. T point is used to sample the Brillouin zone in molecular dynamics 
simulations, because EOS can only be modified within 5% for the selection of higher number 
of k points. 64 helium atoms are included in the cubic supercell. The densities selected in 
our simulations range from 0.1233 to 0.8 g/cm 3 and temperatures between 4 and 50000 
K, which highlight the regime of the principal Hugoniot. All the dynamic simulations are 
lasted for 6000 steps, and the time steps for the integrations of atomic motion are selected 



according to different densities (temperatures) [28]. Then, the subsequent 1000 steps of 



smulation are used to calculate EOS as running averages. 

III. RESULTS AND DISCUSSION 

Different corrections to QMD simulations have already been used to model the ther- 
mophysical properties of WDM, such as the zero point vibrational energy modification for 
hydrogen and deuterium 29|. However, quantitative descriptions of the ionization of atoms 
in the frame of DFT is still lacking, except for the Drude model for aluminium introduced by 
Mazevet et al. 30], but the simple metallic model is not suitable for studying warm dense 
helium. Here, we adopt similar approximations as described in Ref . {3l| , where the effect of 
atomic ionization are considered, to study the EOS of fluid helium. In standard FT-DFT 
molecular dynamic simulations, the ionization energy is excluded, thus the EOS should be 



corrected as 



3l|: 



E = E Q md + Nf3E iom (1) 
P = Pqmd + (1 + P)^, (2) 

where Eqmd arid Pqmd are data obtained from regular QMD calculations, and N is the 
total number of atoms considered in supercells. m# e and ks represent the mass of helium 
atom and Boltzmann constant. The density and temperature are denoted by p and T, 
respectively. /3 stands for the ionization degree, and E ion is the ionization energy. 

In the present work, only the first ionization level is considered, and the ionization degree 
of helium can be evaluated through Saha equation: 

/3 2f2 . Ei on . , . 

W = A^ eXp( -^T ) ' (3) 



4 




FIG. 1: (Color online) Principal Hugoniot up to 130 GPa is shown for the present work (red 
line), where previous theoretical predictions and experimental measurements are also shown for 
comparison. Experiments: Two stage light gas gun data [2] is denoted by open circle; Laser driven 
experimental results 1^] are labelled by open square. Theories: (i) Classical MC simulation results 



321 ] are plotted as green solid line; (ii) ideal plasma's Hugoniot curve [18| is shown as green dotted 



line: (Hi) SCVfl method if- labelled as green dashed line: (iv) ACTEX results 



are blue solid 



line; (v) Previous DFT simulation results jl8l | (with and without accounting for electron excitation 



are shown as black dashed and solid line, respectively; (vi) PIMC results 
dotted line. Inset is the plot of shock wave velocity verse mass velocity. 



181 are labelled as black 



h 2 

X = ^2nm e k B r (4) 
where, m e and Q present the mass of electron and volume of the supercell. The ionization en- 
ergy E ion is determined by E ion = E He+ {p, T) - E He (p, T), where E He+ {p, T) and E He {p, T) 
are total energy of He + and He at the relative density and temperature, respectively. 

Based on the approximations mentioned above, the calculated EOS was examined theo- 
retically along the principal Hugoniot, the locus of states that satisfy the Rankine- Hugoniot 
(RH) equations, which are derived from conservation of mass, momentum, and energy across 
the front of shock waves. The RH equations describe the locus of states in (E, P, \^)-space 
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FIG. 2: Ionization energy and ionization degree are plotted as functions of pressure along the 
principal Hugoniot in the upper and lower panels, respectively. 

satisfying the following relation: 

Ei-E = ^{P 1 + P Q )(V-V ), (5) 

Pi - Po = Pou s u p , (6) 

V l = V (l-u p /u s ), (7) 

where E, P, V present internal energy, pressure, volume, and subscripts and 1 present 
the initial and shocked state, respectively. In Eqs. (jEJ) and (J7J), u s is the velocity of the 
shock wave and u p corresponds to the mass velocity of the material behind the shock front. 
In our present simulations, the initial density for helium is po=0-1233 g/cm 3 and the liquid 
specimen is controlled at a temperature of 4 K, where the internal energy is E =— 0.02 
eV/atom. The initial pressure Po can be treated approximately as zero compared to the 
high pressure of shocked states along the Hugoniot. The Hugoniot points are obtained as 
follows: (i) smooth functions are used to fit the internal energy and pressure in terms of 
temperature at sampled density; (ii) then, Hugoniot points are derived from Eq. (jSJ). 

The principal Hugoniot curve for helium is shown in Fig. [TJ where previous experimental 
measurements and theoretical predictions are also plotted for comparison. For pressures be- 
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FIG. 3: (Color online) Comparison between our calculated shock temperature and previous theo- 
retical predictions. Interatomic potential results are shown as open square [2] and solid line 



DFT 



18] and ACTEX 19] results are plotted as dashed and dotted line, respectively. 



ow 30 GPa, our results show good agreements with the data detected by gas gun experiments 



12l | . At higher pressures, the principal Hugoniot curve shows a maximum compression ratio 
Vmax~6 with the crossover pressure around 100 GPa, as have been reported by laser driven 
experiments [13]. Soon after that, the use of quartz as a shock wave standard by Knudson 
et al. has improve the shock data, and r\ max has been reduced to be 5.1 [14]. Our corrected 
QMD simulation results indicate that, the atomic ionization dominates the characteristic of 
the EOS at P > 30, and the principal Hugoniot shows soften behavior with the increase of 
pressure, then reach its maximum (77 max =5.16) with the corresponding pressure of 106 GPa. 
Then, stiff behavior has been found at P > 110 GPa along the Hugoniot. As has been shown 
in Fig. [21 with the increase of pressure along the Hugoniot, ionization energy decreases, and 
continuous increase of atomic ionization has become considerable in determining the EOS 
of fluid helium. Maximum compression ratio is achieved with the ionization degree of 40%. 
The wide-range behavior of our simulated principal Hugoniot for helium shows excellent 
agreement with experimental ones. 

Concerning other theoretical models, classical MC results show stiff behavior along the 
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Hugoniot, and rj max lies around 4 at 1000 GPa [32]. Direct DFT simulation results show stiff 
behavior up to 110 GPa, while, by combining PIMC simulations, r] max was reported to be 
5.3 at 360 GPa [18j], but the predicted pressure is still too high compared with laser driven 
experiments. ACTEX [19] and SCVH [8] calculations predict rj max to be 6 near 100 GPa 
and 300 GPa respectively, but the results do not agree with those of Knudson et al. jl^| . 

Temperature, which is focused as one of the most important parameters in experiments, is 
difficult to be measured because of the uncertainty in determining the optical-intensity loss 
for ultraviolet part of the spectrum in adiabatic or isentropic shock compressions, especially 



for the temperature exceeding several electron-volt |33j. QMD simulations provide powerful 
tools to predict shock temperature. The calculated Hugoniot temperature has been shown 
in Fig. [3]as a function of pressure along the Hugoniot. Previous theoretical predictions, such 



as interatomic potential 
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15|, DFT-MD (excited electron) |l8|, and ACTEX [19j are also 



provided for comparison. Disagreements have been found to begin at around 30 GPa with 
the relative temperature of 20000 K, which highlights the starting point of ionization on the 
Hugoniot, and the predicted temperatures by those models are higher at a given pressure 
compared with our calculations. 



IV. CONCLUSION 



In summary, We have performed DFT based QMD simulations to study the thermophys- 
ical properties of helium under extreme conditions. The Hugoniot EOS has been evaluated 
through QMD calculations and corrected by taking into account the atomic ionization de- 
scribed by Saha equation, where only first ionization is considered. The corrected Hugoniot 
has been proved accord well with the experimental data in a wide range of shock conditions, 
which thus indicates the importance of atomic ionization. Maximum compression ratio of 
5.16 reveals at round 106 GPa along the principal Hugoniot, and the softened characteristic 
of the Hugoniot has been demonstrated by the contributions from atomic ionization. 
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